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ABSTRACT 

In the presence of strong clustering, low-resolution surveys measure the summed con- 
tributions of groups of sources within the beam. The counts of bright intensity peaks 
are therefore shifted to higher flux levels compared to the counts of individual sources 
detected with high-resolution instruments. If the beam-width corresponds a sizable 
fraction of the clustering size, as in the case of Planck/HFI, one actually detects the 
fluxes of clumps of sources. We argue that the distribution of clump luminosities can 
be modelled in terms of the two- and three-point correlation functions, and apply our 
formalism to the Planck/HFI 850 /im surveys. The effect on counts is found to be large 
and sensitive also to the evolution of the three-point correlation function; in the ex- 
treme case that the latter function is redshift-independent, the source confusion due to 
clustering keeps being important above the canonical 5a detection limit. Detailed sim- 
ulations confirm the reliability of our approach. As the ratio of the beam-width to the 
clustering angular size decreases, the observed fluxes approach those of the brightest 
sources in the beam and the clump formalism no longer applies. However, simulations 
show that also in the case of the Herschel/ SPIRE 500 /xm survey the enhancement of 
the bright source counts due to clustering is important. 
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1 INTRODUCTION 

It has long been known (Eddington 1913) that the counts 
(or flux estimates) of low signal-to-noise sources are biased 
high. Contributions to the noise arise both from the in- 
strument and from source confusion. The latter, which may 
dominate already at relatively bright fluxes in the case of 
low-resolution surveys, comprises the effect of Poisson fluc- 
tuations and of source clustering. For angular scales where 
the correlation of source positions is significant, the ratio 
of clustering to Poisson fluctuations increases with angu- 
lar scale (De Zotti et al. 1996), so that the confusion limit 
can be set by the effect of clustering. This is likely the case 
for some far-IR and for sub-mm surveys from space, due 
to the relatively small primary apertures and to the pres- 
ence of strongly clustered populations (Scott & White 1999; 
Haiman & Knox 1999; Magliocchetti et al. 2001; Perrotta et 
al. 2003; Negrello et al. 2004). 

While the effect of Poisson confusion on source counts 



has been extensively investigated both analytically (Scheuer 
1957; Murdoch et al. 1973; Condon 1974; Hogg & Turner 
1998) and through numerical simulations (Eales et al. 2000; 
Hogg 2001; Blain 2001), the effect of clustering received 
far less attention. The numerical simulations by Hughes & 
Gaztanaga (2000) focused on the sampling variance due to 
clustering. Algorithms successfully simulating the 2D distri- 
bution of clustered sources over sky patches as well as over 
the full sky have been presented by Argiieso et al. (2003) 
and Gonzalez-Nuevo et al. (2004), who also discussed sev- 
eral applications. Some efforts have been made to quantify 
theoretically the effect of clustering on the confusion noise 
but not on the source counts (Toffolatti et al. 1998; Negrello 
et al. 2004; Takeuchi & Ishii 2004). 

In this paper we will use the counts of neighbours for- 
malism (Peebles 1980, hereafter P80) and numerical simu- 
lations to address the effects of clustering on source fluxes 
as measured by low-angular resolution surveys and, conse- 
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quently, on source counts. The formalism is described in 
Section 2, the numerical simulations in Section 3, while in 
Section 4 we present applications to the 850 /jm survey with 
Planck/HFI and to the 500 jun Herschel/SPIRE surveys. In 
Section 5 we summarize our main conclusions. 

We adopt a flat cold dark matter (CDM) cosmology 
with n ,A = 0.7 and h = tfo/lOOkms -1 Mpc -1 = 0.7, 
consistent with the first-year WMAP results (Spergel et al. 
2003). 



2 FORMALISM 

While in the case of a Poisson distribution a source is ob- 
served on top of a background of unresolved sources that 
may be either above or below the all-sky average, in the 
case of clustering sources are preferentially located in over- 
dense environments and one therefore measures the sum of 
all physically related sources in the resolution element of 
the survey (on top of Poisson fluctuations due to unrelated 
sources seen in projection). 

If the beam-width corresponds to a substantial fraction 
of the clustering size, the observed flux is, generally, well 
above that of any individual source. Thus, to estimate the 
counts we cannot refer to the source luminosity function, 
but must define the distribution of luminosities, L c i(z), of 
source "clumps" , as a function of redshift z. It has long been 
known (Kofman et al. 1994; Taylor & Watts 2000; Kayo et 
al. 2001) that a log-normal function is remarkably successful 
in reproducing the statistics of the matter-density distribu- 
tion found in a number of N-body simulations performed 
within the CDM framework, not only in weakly non-linear 
regimes, but also in more strongly non-linear regimes, up 
to density contrasts 8 « 100. Furthermore, it displays the 
correct asymptotic behaviour at very early times, when the 
density field evolves linearly and its distribution is still very 
close to the initial Gaussian one (Coles & Jones 1991). 

If light is a (biased) tracer of mass, fluctuations in the 
luminosity density should obey the same statistics of the 
matter-density field; we therefore adopt a log-normal shape 
for the distribution of L c j. Such a function is completely 
specified by its first (mean) and second (variance) moments, 
that can be evaluated by using the counts of neighbours 
formalism. The mean number of objects inside a volume V 
centered on a given source, {N) p , is [see Eq. (36.23) of P80]: 



(N)p 



[l+£(r)]dV 



(1) 



where n is the mean volume density of the sources and £ 
is their two-point spatial correlation function. The excess of 
objects (with respect to a random distribution) around the 
central one is then given by the second term on the right- 
hand side of this equation. The variance around the mean 
value {N)p can instead be written as [Eq. (36.26) of P80]: 

{(N - (N) p ) 2 ) p = (N} p + 

n f I [C(n, fa) + £(n 2 ) - i{r 1 )i{r 2 )]dV 1 dV 2 . (2) 
Jv Jv 

If the first term on the right-hand side dominates, the vari- 
ance is approximately equal to the mean, as in the case of a 
Poisson distribution. 

The second term on the right-hand side of Eq. (2) is 



related to the skewness of the source distribution and ex- 
hibits a dependence also on the reduced part of the three- 
point angular correlation function, f , for which we adopt the 
standard hierarchical formula: 

C(ri,r a ) = QUr^i^) + Hr^im) + i{r 12 )i{r 2 )]. (3) 

In the local Universe, observational estimates of the am- 
plitude Q indicate nearly constant values, in the range 
Q ~ 0.6-1.3, on scales smaller than ~ 10 Mpc/h (Peebles 
& Groth 1975; Fry & Seldner 1982; Jing & Boerner 1998, 
2004). On larger scales, i.e. in the weakly non-linear regime 
where £ % 1, the non- linear perturbation theory - corrobo- 
rated by N-body simulations - instead predicts Q to show a 
dependence on the scale (Fry 1984; Fry et al. 1993; Jing & 
Boerner 1997; Gaztafiaga & Bernardeau 1998; Scoccimarro 
et al. 1998). 

From Eqs. (1) and (2) we can derive, at a given redshift 
z, the mean luminosity of the clump, L c i, and its variance, 

L c i(z) - L(z) 

+ J dL'K(z)L'$(L',z)- J [l + £(r, z)]dV (4) 
J c Jv 

and 

vl cl (z)=ai+ I dL'K(z)L' 2 $(L',z) [ [1 + £(r, z)}dV 
J c Jv 

+ / dL'K(z)L'$(L',z) ■ 
J c 

[C(ri, ra, z) + £(n 2 , z) - £(n, z)t(ra, z)\dV 1 dV 2 . (5) 



In Eq. (4) L(z) represents the mean luminosity of the sources 
located at redshift z 



= J c dL'K(z)L'<i>(L',z) 
f dL'K(z)<f>(L',z) 



(6) 



and accounts for the fact that we are considering the reso- 
lution elements containing at least one source. The second 
term in the right-hand side of Eq. (4) adds the mean con- 
tribution of neighbours, and approaches zero as the angu- 
lar resolution increases (and correspondingly the volume V 
decreases as the area of the resolution element). In both 
Eq. (4) andEq. (5) K(z) = {l + z)L[v{l + z)]/L{y) is the K- 
correction for monochromatic observations at the frequency 
v and $ is the luminosity function (LF) of the sources. The 
range of integration in luminosity is C =[L m in, £max], where 
Lmin and L max are, respectively, the minimum and the max- 
imum intrinsic luminosity of the sources. The integrals over 
volume are carried out up to r max = 3ro, ro being defined 
by C(ro) = 1. 

The terms on the right-hand side of Eq. (5) account 
for fluctuations around L m (z) (first term), of neighbour lu- 
minosities (second term) and of neighbour numbers (third 
term). The third term has a much steeper dependence on 
the angular resolution than the second one and becomes 
quickly negligible as the area of the resolution element de- 
creases. However, as we will show in Section 4, it is impor- 
tant at the angular resolution of Planck/HFI, implying that 
the distribution of fluxes yielded by such instrument carries 
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information on the evolution of the three-point correlation 
function. 

Suppose that the survey is carried out with an instru- 
ment having a Gaussian angular response function, f(6): 



f(0) = e- (e/e) l\ 

where O relates to the FWHM through: 
FWHM 



e = 



(7) 



(8) 



2\/21n2 

If we adopt the usual power-law model £(r) = (ro/r) 1,8 cut 
at some r max , the integral of the two-point correlation func- 
tion in Eqs. (4) and (5) gives, for a source at redshift z: 



■h 



Z{r)f{6)dV~ 25.9 r^ [D A (z)Qf 



(9) 



where Da(z) is the angular diameter distance and we have 
assumed Da(z)Q < r max . This equation shows that, in the 
case of strongly evolving, highly clustered sources, observed 
with poor angular resolution, the mean number of physically 
correlated neighbours of a galaxy falling within the beam, 
nj?, and therefore their contribution to the observed flux, 
can be quite significant. 

The probability that the sum of luminosities of sources 
in a clump amounts to L c i is then given by: 



p(Lcl,z) 



where 



M«) 



<?g( z ) 



cxp[-l[ln(L ei )- M9 (z)] 2 /q 2 (z)] 
y/2ira1(z)L c i 



In 



In 



Ld(z) 



a 2 L Jz) + L cl (z) 



L c i{z) 



(10) 



(11) 



(12) 



An estimate of the "clump" luminosity function, 
^ clump (Ld,z), is then provided by Eq. (10) apart from the 
normalization factor that we obtain from the conservation 
of the luminosity density: 



*c 



ap (L c i,z)LcidLci = / &(L,z)LdL 



(13) 



The shift to higher luminosities of ^dump compared to $ 
is then compensated by a decrease of the former function 
compared to the latter at low luminosities, as low luminosity 
sources merge to produce a higher luminosity clump. 

If the clustering terms in Eq. (4) and Eq. (5) are much 
larger than those due to individual sources (L m and erf,, 
respectively), the survey will detect "clumps" rather then 
individual sources, and we can use ^cium P (£ci, z) to estimate 
the counts. 

Both theoretical arguments and observational data indi- 
cate that the positions of powerful far-IR galaxies detected 
by SCUBA surveys ("SCUBA galaxies") are highly corre- 
lated (see e.g. Smail et al. 2003; Negrello et al. 2004; Blain et 
al. 2004) so that their confusion fluctuations arc dominated 
by clustering effects. On the contrary, Poisson fluctuations 
dominate in the case of the other extragalactic source pop- 
ulations contributing to the sub-millimeter counts (spiral 
and starburst galaxies, whose clustering is relatively weak; 



cf. e.g. Madgwick et al. 2003). Therefore we will neglect the 
clustering of the latter populations (which are however in- 
cluded in the estimates of source counts) and will apply the 
above formalism to SCUBA galaxies only; for their evolv- 
ing two-point spatial correlation function, £(r, 2), we adopt 
model 2 of Negrello et al. (2004). According to such model 



£(r,z)=b (M eS ,z)£nM(r,z), 



(14) 



where £,dm is the non-linear two-point spatial correlation 
function of dark matter, computed with the recipe by Pea- 
cock & Dodds (1996; see also Smith et al. 2003), adopting 
a CDM spectrum for the primordial density perturbations, 
with an index n = 1, a shape parameter V = 0.2 and a nor- 
malization erg = 0.8 (see e.g. Lahav et al. 2002; Spergel et al. 
2003); b(M e s,z) is the redshift-dependent (linear) bias fac- 
tor (Sheth & Tormen 1999), M c g being the effective mass of 
the dark matter halos in which SCUBA sources reside. Fol- 
lowing Negrello et al. (2004) we set M cB = 1.8 x 10 13 M©/h, 
which yields values of £(r, z) consistent with the available 
observational estimates. Note that a one-to-one correspon- 
dence between haloes and sources has been assumed. 

For the amplitude, Q, of the three-point angular corre- 
lation function, £ [Eq. (3)], we consider three models: 

(i) Q(z) = Q(0) = 1, 

(ii) Q{z) = Q{0)/b(M eB ,z) = l/b(M cB ,z), 
(iii) Q{z) = Q(0)/b 2 (M cS ,z) = l/b 2 (M cB ,z), 

where we have neglected any dependence of Q on scale. Both 
calculations based on perturbation theory (e.g. Juszkiewicz, 
Bouchet & Colombi 1993; Bernardeau 1994) and N-body 
simulations (e.g. Colombi, Bouchet & Hernquist 1996; Sza- 
pudi et al. 1996) suggest that model (i) applies to dark 
matter. On the other hand, the three-point correlation of 
luminous objects decreases as the bias factor b increases 
(Bernardeau & Schaeffer 1992, 1999; Szapudi et al. 2001); 
the formula (ii), derived from perturbation theory, is ex- 
pected to hold on scales > 10 Mpc/h (Fry & Gaztafiaga 
1993) while for scales smaller than these, Szapudi et al. 
(2001) quote model (iii) as a phenomenological rule derived 
from N-body simulations. A more realistic model should al- 
low for the dependence of Q on the linear scale, induced by 
the increasing strength of non-linear effects with decreas- 
ing scale. The models (i)-(iii) may thus bracket the true 
behaviour of Q(z), which should, however, be not far from 
model (iii) for the scales of interest here. 



3 NUMERICAL SIMULATIONS 

To test the reliability of our analytical approach, numerical 
simulations of sky patches were performed using the fast al- 
gorithm recently developed by Gonzalez- Nuevo et al. (2004) . 
Only SCUBA galaxies were taken into account. Sources were 
first randomly distributed over the patch area, with surface 
densities given, as a function of the flux density, by the model 
of Granato et al. (2004); we have considered sources down to 
Smin = 0.01 mjy. Then, the projected density contrast as a 
function of position, <5(x), was derived and its Fourier trans- 
form, <5(k), was computed. Next, in Fourier space, we added 
the power spectrum corresponding to the two-point angular 
correlation function, w(9) , given by model 2 of Negrello et al. 
(2004), to the white noise power spectrum corresponding to 
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Figure 1. Differential source counts of SCUBA galaxies at 
850 /im (solid line) compared with counts expected in the 
case of observations performed with the angular resolution of 
Planck/HFI for the three models for the evolution of the am- 
plitude of the three-point correlation function, Q(z), (model (i): 
short dashes; model (ii): long dashes; model (iii): dot-dashed line). 
The summed contributions from (unclustered) spiral galaxies, 
starburst galaxies, and extragalactic radio sources is represented 
by the dotted line. The filled circles with error bars show the 
counts estimated from the simulations (see text). The vertical dot- 
ted line shows the Planck/HFI 5cr detection limit (Sd = 288 mJy) 
estimated by Negrello et al. (2004) allowing for clustering fluctu- 
ations. 



the initial spatial distribution and obtained the transformed 
density field of spatially correlated sources. Then we apply 
the inverse Fourier transform to get the projected distribu- 
tion of clustered sources in the real space and we associate 
randomly the fluxes to the simulated sources according to 
the differential counts predicted by the model of Granato et 
al. (2004) (for more details, see Gonzalez-Nuevo et al. 2004) 
Simulations were carried out for surveys in: 

> the 850 /im channel of the High Frequency Instrument 
(HFI) of the ESA Planck satellite (FWHM=300"; Lamarre 
et al. 2003); 

> the 500 /im channel of the Spectral and Photometric 
Imaging REceiver (SPIRE) of the ESA Herschel satellite 
(FWHM=34.6"; Griffin et al. 2000). 

In the Planck/HFI case, we have simulated sky patches of 
12°. 8 x 12°.8deg 2 with pixels of 1.5 x 1.5arcmin 2 . In the 
Herschel/SPIRE case the patches were of 3°. 28 x 3°.28deg 2 
and the pixel size was 1/3 of the FWHM. 

To check the simulation procedure we have compared 
the confusion noise, o, obtained from them with the an- 
alytical results by Negrello et al. (2004). In the case of 
Planck/HFI, from the simulations we get op — 20mJy 
and ac = 50mJy for Poisson distributed and clustered 
sources, respectively, in very good agreement with the ana- 
lytic results, op — 20 mJy and oc — 53 mJy. For the Her- 
schel/SPIRE channel the simulations give op — 6.4 mJy 
and oc = 1.2 mJy, to be compared with op — 6.0 mJy and 
oc = 3.2 mJy obtained by Negrello et al. (2004). The appar- 
ent discrepancy of the results for oc actually corresponds to 



the uncertainty in this quantity, very difficult to determine 
from the simulations when Poisson fluctuations dominate. 

Note that for the very steep counts predicted by the 
Granato et al. (2004) model, that accurately matches the 
SCUBA and MAMBO data, both the Poisson and the clus- 
tering fluctuations, for the angular resolutions considered 
here, are independent of the flux limit. We have checked 
that the results obtained from simulations are independent 
of the pixel size used. 



4 RESULTS AND DISCUSSION 

The results for the Planck/HFI 850 /^m channel are shown in 
Fig. 1, where the solid line gives the counts of SCUBA galax- 
ies predicted by the physically grounded model of Granato 
et al. (2004), and the dotted line gives the summed counts 
of spiral and starburst galaxies, and of extragalactic radio 
sources. For the redshift-dependent luminosity functions of 
spiral and starburst galaxies we have adopted the same phe- 
nomenological models as in Negrello et al. (2004), while for 
radio galaxies we have used the model by De Zotti et al. 
(2004). The short-dashed, long-dashed, and dot-dashed lines 
show the counts of "clumps" expected from our analytic 
formalism for the three evolution models of the amplitude 
Q of the three-point correlation function. As expected, at 
the Planck/HFI resolution the counts of "clumps" are sensi- 
tive to the evolution of the three-point correlation function. 
Clearly the predicted counts below ~ 50 mJy, where Poisson 
fluctuations become important, are of no practical use; they 
are shown just to illustrate how the formalism accounts for 
the disappearance of lower luminosity objects which merge 
into the "clumps". 

The filled circles with error bars in Fig. 1 are obtained 
filtering the simulated maps with a Gaussian response func- 
tion of 5' FWHM to mimic Planck/HFI observations. The 
lower flux limit of the estimated counts is set by Poisson 
fluctuations. It should be noted that the procedure used for 
the simulations takes into account only the two-point angu- 
lar correlation function and does not allow us to deal with 
the three-point correlation function and with its cosmolog- 
ical evolution, which are included in the analytic model. In 
principle it is possible to go the other way round, i.e. to eval- 
uate from the simulations the reduced angular bispectrum, 
b r , by applying the standard Fourier analysis (Gonzalez- 
Nuevo et al. 2004; Argiieso et al. 2003) and infer from it 
an estimate of the three-point correlation function weighted 
over the redshift distribution. However, the relationship of 
the three-point correlation function with the angular bispec- 
trum is through a six dimensional integral which is really 
difficult to deal with in practice (Szapudi 2004) . An alterna- 
tive possibility consists in computing the number of triplets 
(above a fixed flux threshold) and in using the estimator de- 
veloped by Szapudi & Szalay (1998) to derive an estimate of 
the three-point angular correlation function. On the whole, 
these methods turn out to be impractical to take into ac- 
count also the effect of the evolving three-point correlation 
function when comparing simulations with analytic results. 
On the other hand, it is very reassuring that the counts ob- 
tained from simulations are within the range spanned by 
analytic models. 

In the Herschel case, the beam encompasses only a small 
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Figure 2. Left-hand panel: same as in Fig. 1 but for the 500 (im channel of Herschel/SPIRE. The results for models (ii) and (iii) for 
the evolution of Q(z) overlap, while model (i) is only slightly higher. Bight-hand panel: ratio between the differential counts of SCUBA 
galaxies estimated from the simulations for the Herschel/SPIRE resolution and those predicted by the model of Granato et al. (2004). The 
dotted vertical line, in both panels, corresponds to the 5<r detection limit (S^ = 39 mJy) estimated by Negrello et al. (2004) accounting 
also for clustering fluctuations. 



fraction of the "clump" and therefore the observed flux is 
generally dominated by the single brightest source in the 
beam. Thus, the analytic model described in Section 2 is no 
longer applicable. As illustrated by the left-hand panel of 
Fig. 2, such formalism would strongly over-predict the ob- 
served counts at bright flux densities (and the results are 
essentially independent of the three-point correlation func- 
tion). On the other hand, the simulations (filled circles with 
error bars) show that neighbour sources appreciably con- 
tribute to the observed fluxes, hence to the counts at bright 
flux density levels, as more clearly illustrated by the right- 
hand panel of Fig. 2. Such count estimates were obtained 
filtering the simulated maps with a Gaussian response func- 
tion of FWHM=34.6", appropriate for the Herschel 500 ^m 
channel; again, the lower flux limit of the estimated counts 
is set by Poisson fluctuations. 



5 CONCLUSIONS 

Theoretical arguments and observational data converge in 
indicating that the very luminous (sub)-mm sources de- 
tected by SCUBA and MAMBO surveys are highly clustered 
(clustering radius ro ~ Sh" 1 Mpc). On the other hand, the 
limited sizes of telescopes of forthcoming space instruments 
operating in the sub-mm domain, such as Planck/HFI and 
Herschel/SPIRE, imply relatively poor angular resolutions 
(FWHM ~ 5' for Planck at 850 fira, and ~ 34.6" for Her- 
schel at at 500 ^m) . 

In the Planck/HFI case the summed fluxes of physically 
correlated sources within the beam are generally higher than 
the luminosity of the brightest source in the beam, so that 
the outcome of the surveys will be counts of "clumps" of 
sources rather than of individual sources. We have argued 
that the luminosity distribution of such "clumps" at any red- 



shift can be modelled with a log-normal function, with mean 
determined by the average source luminosity and by the av- 
erage of summed luminosities of neighbours, and variance 
made of three contributions (the variances of source lumi- 
nosities, of neighbour luminosities, and of neighbour num- 
bers). The latter contribution depends on the three-point 
correlation function, so that the counts of "clumps" provide 
information on this elusive quantity and on its cosmologi- 
cal evolution. Under the, rather extreme, assumption that 
the coefficient, Q, of the three-point correlation function is 
independent of redshift, the counts of "clumps" extend be- 
yond the formal 5<r detection limit. In the more likely cases 
of Q decreasing as b^ 1 or 6~ 2 , b being the bias factor, the 
"clumps" will only show up in Planck maps as < 5<j fluctu- 
ations. Anyway, the Planck surveys will provide a rich cat- 
alogue of candidate proto-clusters at substantial redshifts 
(typically at z ~ 2-3), very important to investigate the 
formation of large scale structure and, particularly, to con- 
strain the evolution of the dark energy thought to control 
the dynamics of the present day universe. Detailed numeri- 
cal simulations carried out using the fast algorithm recently 
developed by Gonzalez-Nuevo et al. (2004) are fully consis- 
tent with the analytic results, although a full comparison 
would require an upgrade of the algorithm to include the 
effect of the evolving three-point correlation function. 

As the ratio of the beam-width to the clustering an- 
gular size decreases, the observed fluxes approach those of 
the brightest sources in the beam and the "clump" formal- 
ism no longer applies. However, simulations show that also 
in the case of the Herschel/SPIRE 500 /im survey the con- 
tribution of neighbours to the observed fluxes enhances the 
bright tail of the observed counts. Due to the extreme steep- 
ness of such tail, as predicted by the model of Granato et al. 
(2004), even a modest addition to the fluxes of the bright- 
est sources may lead to counts at flux densities ~ 100 mJy 
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several times higher than would be observed with a high 
resolution instrument. It should be noted that, in the case 
of strong clustering, the canonical 5<r detection limit (shown 
by the vertical dotted line in both panels of Fig. 2) does not 
frees the observed counts from the confusion bias. 
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ABSTRACT 

In the presence of strong clustering, low-resolution surveys measure the summed con- 
tributions of groups of sources within the beam. The counts of bright intensity peaks 
are therefore shifted to higher flux levels compared to the counts of individual sources 
detected with high- resolution instruments. If the beam- width corresponds a sizable 
fraction of the clustering size, as in the case of Planck/HFI, one actually detects the 
fluxes of clumps of sources. We argue that the distribution of clump luminosities can 
be modelled in terms of the two- and three-point correlation functions, and apply our 
formalism to the Planck/HFI 850 /tm surveys. The effect on counts is found to be large 
and sensitive also to the evolution of the three-point correlation function; in the ex- 
treme case that the latter function is redshift-independent, the source confusion due to 
clustering keeps being important above the canonical 5<r detection limit. Detailed sim- 
ulations confirm the reliability of our approach. As the ratio of the beam-width to the 
clustering angular size decreases, the observed fluxes approach those of the brightest 
sources in the beam and the clump formalism no longer applies. However, simulations 
show that also in the case of the Herschel/SPIRE 500 /tm survey the enhancement of 
the bright source counts due to clustering is important. 

Key words: galaxies: evolution, counts - sub-millimetre - clustering: models 



1 INTRODUCTION 

It has long been known (Eddington 1913) that the counts 
(or flux estimates) of low signal-to-noise sources are biased 
high. Contributions to the noise arise both from the in- 
strument and from source confusion. The latter, which may 
dominate already at relatively bright fluxes in the case of 
low-resolution surveys, comprises the effect of Poisson fluc- 
tuations and of source clustering. For angular scales where 
the correlation of source positions is significant, the ratio 
of clustering to Poisson fluctuations increases with angu- 
lar scale (De Zotti et al. 1996), so that the confusion limit 
can be set by the effect of clustering. This is likely the case 
for some far-IR and for sub-mm surveys from space, due 
to the relatively small primary apertures and to the pres- 
ence of strongly clustered populations (Scott & White 1999; 
Haiman & Knox 1999; Magliocchetti et al. 2001; Perrotta et 
al. 2003; Negrello et al. 2004). 

While the effect of Poisson confusion on source counts 



has been extensively investigated both analytically (Scheuer 
1957; Murdoch et al. 1973; Condon 1974; Hogg & Turner 
1998) and through numerical simulations (Eales et al. 2000; 
Hogg 2001; Blain 2001), the effect of clustering received 
far less attention. The numerical simulations by Hughes & 
Gaztafiaga (2000) focused on the sampling variance due to 
clustering. Algorithms successfully simulating the 2D distri- 
bution of clustered sources over sky patches as well as over 
the full sky have been presented by Argiieso et al. (2003) 
and Gonzalez-Nuevo et al. (2004), who also discussed sev- 
eral applications. Some efforts have been made to quantify 
theoretically the effect of clustering on the confusion noise 
but not on the source counts (Toffolatti et al. 1998; Negrello 
et al. 2004; Takeuchi & Ishii 2004). 

In this paper we will use the counts of neighbours for- 
malism (Peebles 1980, hereafter P80) and numerical simu- 
lations to address the effects of clustering on source fluxes 
as measured by low-angular resolution surveys and, conse- 



© 0000 RAS 



quently, on source counts. The formalism is described in 
Section 2, the numerical simulations in Section 3, while in 
Section 4 we present applications to the 850 fim survey with 
Planck/HFI and to the 500 pm Herschel/SPIRE surveys. In 
Section 5 we summarize our main conclusions. 

We adopt a flat cold dark matter (CDM) cosmology 
with O ,A = 0.7 and h = ffo/lOOkms" 1 Mpc" 1 = 0.7, 
consistent with the first-year WMAP results (Spergel et al. 
2003). 



2 FORMALISM 

While in the case of a Poisson distribution a source is ob- 
served on top of a background of unresolved sources that 
may be either above or below the all-sky average, in the 
case of clustering sources are preferentially located in over- 
dense environments and one therefore measures the sum of 
all physically related sources in the resolution element of 
the survey (on top of Poisson fluctuations due to unrelated 
sources seen in projection). 

If the beam-width corresponds to a substantial fraction 
of the clustering size, the observed flux is, generally, well 
above that of any individual source. Thus, to estimate the 
counts we cannot refer to the source luminosity function, 
but must define the distribution of luminosities, L c i(z), of 
source "clumps" , as a function of redshift z. It has long been 
known (Kofman et al. 1994; Taylor & Watts 2000; Kayo et 
al. 2001) that a log-normal function is remarkably successful 
in reproducing the statistics of the matter-density distribu- 
tion found in a number of N-body simulations performed 
within the CDM framework, not only in weakly non-linear 
regimes, but also in more strongly non-linear regimes, up 
to density contrasts S « 100. Furthermore, it displays the 
correct asymptotic behaviour at very early times, when the 
density field evolves linearly and its distribution is still very 
close to the initial Gaussian one (Coles & Jones 1991). 

If light is a (biased) tracer of mass, fluctuations in the 
luminosity density should obey the same statistics of the 
matter-density field; we therefore adopt a log-normal shape 
for the distribution of L c i. Such a function is completely 
specified by its first (mean) and second (variance) moments, 
that can be evaluated by using the counts of neighbours 
formalism. The mean number of objects inside a volume V 
centered on a given source, {N) p , is [see Eq. (36.23) of P80]: 



(N) P 



n I [l+£(r)]dV 
'v 



(1) 



where n is the mean volume density of the sources and £ 
is their two-point spatial correlation function. The excess of 
objects (with respect to a random distribution) around the 
central one is then given by the second term on the right- 
hand side of this equation. The variance around the mean 
value (N)p can instead be written as [Eq. (36.26) of P80]: 

{(N - {N) p f) p = {N) p + 

n / [C(ri,r 2 )+$(ri2)-$(ri)$(r 2 )]dVidV2 ■ (2) 

Jv Jv 

If the first term on the right-hand side dominates, the vari- 
ance is approximately equal to the mean, as in the case of a 
Poisson distribution. 

The second term on the right-hand side of Eq. (2) is 



related to the skewness of the source distribution and ex- 
hibits a dependence also on the reduced part of the three- 
point angular correlation function, £, for which we adopt the 
standard hierarchical formula: 



C(n, r 2 ) = QK(ri)£(r 2 ) + £(n)£(ri 2 ) + £(r 12 )£(r 2 )]. 



(3) 



In the local Universe, observational estimates of the am- 
plitude Q indicate nearly constant values, in the range 
Q ~ 0.6-1.3, on scales smaller than ~ 10 Mpc/h (Peebles 
& Groth 1975; Fry & Seldner 1982; Jing & Boerner 1998, 
2004). On larger scales, i.e. in the weakly non-linear regime 
where £ £ 1, the non-linear perturbation theory - corrobo- 
rated by N-body simulations - instead predicts Q to show a 
dependence on the scale (Fry 1984; Fry et al. 1993; Jing & 
Boerner 1997; Gaztafiaga & Bernardeau 1998; Scoccimarro 
et al. 1998). 

From Eqs. (1) and (2) we can derive, at a given redshift 
z, the mean luminosity of the clump, L c i, and its variance, 



< 



L c i(z) = L(z) 



+ dL'K(z)L'$(L',z)- [l+Z(r,z)]dV (4) 

J £ JV 



and 



<rl cl {z) = al+ I dL'K(z)L' 2 $(L',z) / [l + £(r,z)]dV 
Jc Jv 

2 



+ 



/ dL'K{z)L'${L',z) 



I f [ar u r 2 ,z)+i{r l2 ,z)-i{r u z)i{r 2 ,z)]dV,dV 2 . (5) 
Jv Jv 

In Eq. (4) L(z) represents the mean luminosity of the sources 
located at redshift z 



L(z) 



f c dL'K(z)L'*(L',z) 

J £ dL>K(z)$(L',z) 



(6) 



and accounts for the fact that we are considering the reso- 
lution elements containing at least one source. The second 
term in the right-hand side of Eq. (4) adds the mean con- 
tribution of neighbours, and approaches zero as the angu- 
lar resolution increases (and correspondingly the volume V 
decreases as the area of the resolution element). In both 
Eq. (4) andEq. (5) K(z) = (l + z)L[v(l + z)]/L(v) is the K- 
correction for monochromatic observations at the frequency 
v and $ is the luminosity function (LF) of the sources. The 
range of integration in luminosity is L =[L m in, imm], where 
Z/min and Z/max are, respectively, the minimum and the max- 
imum intrinsic luminosity of the sources. The integrals over 
volume are carried out up to r max = 3ro, ro being defined 
by £(r„) = 1. 

The terms on the right-hand side of Eq. (5) account 
for fluctuations around L m (z) (first term), of neighbour lu- 
minosities (second term) and of neighbour numbers (third 
term). The third term has a much steeper dependence on 
the angular resolution than the second one and becomes 
quickly negligible as the area of the resolution element de- 
creases. However, as we will show in Section 4, it is impor- 
tant at the angular resolution of Planck/HFI, implying that 
the distribution of fluxes yielded by such instrument carries 
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information on the evolution of the three-point correlation 
function. 

Suppose that the survey is carried out with an instru- 
ment having a Gaussian angular response function, /(#): 



f(0) 



-(«/e) 2 /2 



where G relates to the FWHM through: 



e 



FWHM 



2%/21n2 



(7) 



(8) 



If we adopt the usual power-law model £(r) = (ro/r) 1 ' 8 cut 
at some r max , the integral of the two-point correlation func- 
tion in Eqs. (4) and (5) gives, for a source at redshift z: 



I 



J2 = / t(r)f{8)dV ~ 25.9 rj 8 [D A (z)e] 1 



(9) 



where Da(z) is the angular diameter distance and we have 
assumed Da{z)Q < r max . This equation shows that, in the 
case of strongly evolving, highly clustered sources, observed 
with poor angular resolution, the mean number of physically 
correlated neighbours of a galaxy falling within the beam, 
nJ2, and therefore their contribution to the observed flux, 
can be quite significant. 

The probability that the sum of luminosities of sources 
in a clump amounts to L c i is then given by: 



p(L c i,z) 
where 



CT s(«) 



exp[-±[HL cl ) - ^(z)} 2 /a 2 g (z)] 
y / 2ira1(z)L c i 



In 



In 



L c i(z) 



a\Jz) + L' cl {z) 
+ 1 



<■(*) 



L c i(z 



(10) 



(11) 



(12) 



An estimate of the "clump" luminosity function, 
*ciump(ioi, z), is then provided by Eq. (10) apart from the 
normalization factor that we obtain from the conservation 
of the luminosity density: 



h 



ap(L c i, z)L c [dL c 



/ $(L,z)LdL 



(13) 



The shift to higher luminosities of * c ium P compared to $ 
is then compensated by a decrease of the former function 
compared to the latter at low luminosities, as low luminosity 
sources merge to produce a higher luminosity clump. 

If the clustering terms in Eq. (4) and Eq. (5) are much 
larger than those due to individual sources (L m and a\, 
respectively), the survey will detect "clumps" rather then 
individual sources, and we can use \]/ c i um p(L c (, z) to estimate 
the counts. 

Both theoretical arguments and observational data indi- 
cate that the positions of powerful far-IR galaxies detected 
by SCUBA surveys ("SCUBA galaxies") are highly corre- 
lated (see e.g. Smail et al. 2003; Negrello et al. 2004; Blain et 
al. 2004) so that their confusion fluctuations are dominated 
by clustering effects. On the contrary, Poisson fluctuations 
dominate in the case of the other extragalactic source pop- 
ulations contributing to the sub-millimeter counts (spiral 
and starburst galaxies, whose clustering is relatively weak; 



cf. e.g. Madgwick et al. 2003). Therefore we will neglect the 
clustering of the latter populations (which are however in- 
cluded in the estimates of source counts) and will apply the 
above formalism to SCUBA galaxies only; for their evolv- 
ing two-point spatial correlation function, £(r, z), we adopt 
model 2 of Negrello et al. (2004). According to such model 



£(r, z) = b (M eS ,z)^oM(r,z), 



(14) 



where £dm is the non-linear two-point spatial correlation 
function of dark matter, computed with the recipe by Pea- 
cock & Dodds (1996; see also Smith et al. 2003), adopting 
a CDM spectrum for the primordial density perturbations, 
with an index n = 1, a shape parameter T = 0.2 and a nor- 
malization <r 8 = 0.8 (see e.g. Lahav et al. 2002; Spergel et al. 
2003); b(M e g,z) is the redshift-dependent (linear) bias fac- 
tor (Sheth & Tormen 1999), M e s being the effective mass of 
the dark matter halos in which SCUBA sources reside. Fol- 
lowing Negrello et al. (2004) we set M eff = 1.8 x 10 13 M /h, 
which yields values of £(r, z) consistent with the available 
observational estimates. Note that a one-to-one correspon- 
dence between haloes and sources has been assumed. 

For the amplitude, Q, of the three-point angular corre- 
lation function, £ [Eq. (3)], we consider three models: 



(i) Q(z) 
(ii) Q(z) 



0(0) = i, 

Q(0)/b(M eS ,z) 



(iii) Q(z) = Q(0)/b 2 (M eS ,z) 



l/b(M eB ,z), 
= l/b 2 (M eB ,z) 



where we have neglected any dependence of Q on scale. Both 
calculations based on perturbation theory (e.g. Juszkiewicz, 
Bouchet & Colombi 1993; Bernardeau 1994) and N-body 
simulations (e.g. Colombi, Bouchet & Hernquist 1996; Sza- 
pudi et al. 1996) suggest that model (i) applies to dark 
matter. On the other hand, the three-point correlation of 
luminous objects decreases as the bias factor b increases 
(Bernardeau & Schaeffer 1992, 1999; Szapudi et al. 2001); 
the formula (ii), derived from perturbation theory, is ex- 
pected to hold on scales > 10 Mpc/h (Fry & Gaztanaga 
1993) while for scales smaller than these, Szapudi et al. 
(2001) quote model (iii) as a phenomenological rule derived 
from N-body simulations. A more realistic model should al- 
low for the dependence of Q on the linear scale, induced by 
the increasing strength of non-linear effects with decreas- 
ing scale. The models (i)-(iii) may thus bracket the true 
behaviour of Q{z), which should, however, be not far from 
model (iii) for the scales of interest here. 



3 NUMERICAL SIMULATIONS 

To test the reliability of our analytical approach, numerical 
simulations of sky patches were performed using the fast al- 
gorithm recently developed by Gonzalez-Nuevo et al. (2004). 
Only SCUBA galaxies were taken into account. Sources were 
first randomly distributed over the patch area, with surface 
densities given, as a function of the flux density, by the model 
of Granato et al. (2004); we have considered sources down to 
£ m in = 0.01 mjy. Then, the projected density contrast as a 
function of position, <5(x), was derived and its Fourier trans- 
form, <5(k), was computed. Next, in Fourier space, we added 
the power spectrum corresponding to the two-point angular 
correlation function, w(0), given by model 2 of Negrello et al. 
(2004), to the white noise power spectrum corresponding to 
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Figure 1. Differential source counts of SCUBA galaxies at 
850 jum (solid line) compared with counts expected in the 
case of observations performed with the angular resolution of 
Planck/HFI for the three models for the evolution of the am- 
plitude of the three-point correlation function, Q(z), (model (i): 
short dashes; model (ii): long dashes; model (iii): dot-dashed line). 
The summed contributions from (unclustered) spiral galaxies, 
starburst galaxies, and extragalactic radio sources is represented 
by the dotted line. The filled circles with error bars show the 
counts estimated from the simulations (see text). The vertical dot- 
ted line shows the Planck/HFI 5<r detection limit (S4 = 288 mJy) 
estimated by Negrello et al. (2004) allowing for clustering fluctu- 
ations. 



the initial spatial distribution and obtained the transformed 
density field of spatially correlated sources. Then we apply 
the inverse Fourier transform to get the projected distribu- 
tion of clustered sources in the real space and we associate 
randomly the fluxes to the simulated sources according to 
the differential counts predicted by the model of Granato et 
al. (2004) (for more details, see Gonzalez-Nuevo et al. 2004) 
Simulations were carried out for surveys in: 

> the 850 |Um channel of the High Frequency Instrument 
(HFI) of the ESA Planck satellite (FWHM=300"; Lamarre 
et al. 2003); 

> the 500 /um channel of the Spectral and Photometric 
Imaging REceiver (SPIRE) of the ESA Herschel satellite 
(FWHM=34.6"; Griffin et al. 2000). 

In the Planck/HFI case, we have simulated sky patches of 
12°. 8 x 12°.8deg 2 with pixels of 1.5 x 1.5arcmin 2 . In the 
Herschel/SPIRE case the patches were of 3°. 28 x 3°. 28 deg 2 
and the pixel size was 1/3 of the FWHM. 

To check the simulation procedure we have compared 
the confusion noise, a, obtained from them with the an- 
alytical results by Negrello et al. (2004). In the case of 
Planck/HFI, from the simulations we get ap = 20 mJy 
and ac = 50 mJy for Poisson distributed and clustered 
sources, respectively, in very good agreement with the ana- 
lytic results, ap = 20 mJy and ac = 53 mJy. For the Her- 
schel/SPIRE channel the simulations give ap = 6.4 mJy 
and ac = 1.2 mJy, to be compared with ap = 6.0 mJy and 
ac = 3.2 mJy obtained by Negrello et al. (2004). The appar- 
ent discrepancy of the results for ac actually corresponds to 



the uncertainty in this quantity, very difficult to determine 
from the simulations when Poisson fluctuations dominate. 

Note that for the very steep counts predicted by the 
Granato et al. (2004) model, that accurately matches the 
SCUBA and MAMBO data, both the Poisson and the clus- 
tering fluctuations, for the angular resolutions considered 
here, are independent of the flux limit. We have checked 
that the results obtained from simulations are independent 
of the pixel size used. 



4 RESULTS AND DISCUSSION 

The results for the Planck/HFI 850 pm channel are shown in 
Fig. 1, where the solid line gives the counts of SCUBA galax- 
ies predicted by the physically grounded model of Granato 
et al. (2004), and the dotted line gives the summed counts 
of spiral and starburst galaxies, and of extragalactic radio 
sources. For the redshift-dependent luminosity functions of 
spiral and starburst galaxies we have adopted the same phe- 
nomenological models as in Negrello et al. (2004), while for 
radio galaxies we have used the model by De Zotti et al. 
(2004). The short-dashed, long-dashed, and dot-dashed lines 
show the counts of "clumps" expected from our analytic 
formalism for the three evolution models of the amplitude 
Q of the three-point correlation function. As expected, at 
the Planck/HFI resolution the counts of "clumps" are sensi- 
tive to the evolution of the three-point correlation function. 
Clearly the predicted counts below ~ 50 mJy, where Poisson 
fluctuations become important, are of no practical use; they 
are shown just to illustrate how the formalism accounts for 
the disappearance of lower luminosity objects which merge 
into the "clumps" . 

The filled circles with error bars in Fig. 1 are obtained 
filtering the simulated maps with a Gaussian response func- 
tion of 5' FWHM to mimic Planck/HFI observations. The 
lower flux limit of the estimated counts is set by Poisson 
fluctuations. It should be noted that the procedure used for 
the simulations takes into account only the two-point angu- 
lar correlation function and does not allow us to deal with 
the three-point correlation function and with its cosmolog- 
ical evolution, which are included in the analytic model. In 
principle it is possible to go the other way round, i.e. to eval- 
uate from the simulations the reduced angular bispectrum, 
b r , by applying the standard Fourier analysis (Gonzalez- 
Nuevo et al. 2004; Argiieso et al. 2003) and infer from it 
an estimate of the three-point correlation function weighted 
over the redshift distribution. However, the relationship of 
the three-point correlation function with the angular bispec- 
trum is through a six dimensional integral which is really 
difficult to deal with in practice (Szapudi 2004). An alterna- 
tive possibility consists in computing the number of triplets 
(above a fixed flux threshold) and in using the estimator de- 
veloped by Szapudi & Szalay (1998) to derive an estimate of 
the three-point angular correlation function. On the whole, 
these methods turn out to be impractical to take into ac- 
count also the effect of the evolving three-point correlation 
function when comparing simulations with analytic results. 
On the other hand, it is very reassuring that the counts ob- 
tained from simulations are within the range spanned by 
analytic models. 

In the Herschel case, the beam encompasses only a small 
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Figure 2. Left-hand panel: same as in Fig. 1 but for the 500 /im channel of Herschel/SPIRE. The results for models (ii) and (iii) for 
the evolution of Q(z) overlap, while model (i) is only slightly higher. Right-hand panel: ratio between the differential counts of SCUBA 
galaxies estimated from the simulations for the Herschel/SPIRE resolution and those predicted by the model of Granato et al. (2004). The 
dotted vertical line, in both panels, corresponds to the 5<r detection limit (S4 = 39mjy) estimated by Negrello et al. (2004) accounting 
also for clustering fluctuations. 



fraction of the "clump" and therefore the observed flux is 
generally dominated by the single brightest source in the 
beam. Thus, the analytic model described in Section 2 is no 
longer applicable. As illustrated by the left-hand panel of 
Fig. 2, such formalism would strongly over-predict the ob- 
served counts at bright flux densities (and the results are 
essentially independent of the three-point correlation func- 
tion). On the other hand, the simulations (filled circles with 
error bars) show that neighbour sources appreciably con- 
tribute to the observed fluxes, hence to the counts at bright 
flux density levels, as more clearly illustrated by the right- 
hand panel of Fig. 2. Such count estimates were obtained 
filtering the simulated maps with a Gaussian response func- 
tion of FWHM=34.6", appropriate for the Herschel 500 pm 
channel; again, the lower flux limit of the estimated counts 
is set by Poisson fluctuations. 



5 CONCLUSIONS 

Theoretical arguments and observational data converge in 
indicating that the very luminous (sub)-mm sources de- 
tected by SCUBA and MAMBO surveys are highly clustered 
(clustering radius ro ~ 8h _1 Mpc). On the other hand, the 
limited sizes of telescopes of forthcoming space instruments 
operating in the sub-mm domain, such as Planck/HFI and 
Herschel/SPIRE, imply relatively poor angular resolutions 
(FWHM ~ 5' for Planck at 850 (im, and ~ 34.6" for Her- 
schel at at 500 jim). 

In the Planck/HFI case the summed fluxes of physically 
correlated sources within the beam are generally higher than 
the luminosity of the brightest source in the beam, so that 
the outcome of the surveys will be counts of "clumps" of 
sources rather than of individual sources. We have argued 
that the luminosity distribution of such "clumps" at any red- 



shift can be modelled with a log-normal function, with mean 
determined by the average source luminosity and by the av- 
erage of summed luminosities of neighbours, and variance 
made of three contributions (the variances of source lumi- 
nosities, of neighbour luminosities, and of neighbour num- 
bers). The latter contribution depends on the three-point 
correlation function, so that the counts of "clumps" provide 
information on this elusive quantity and on its cosmologi- 
cal evolution. Under the, rather extreme, assumption that 
the coefficient, Q, of the three-point correlation function is 
independent of redshift, the counts of "clumps" extend be- 
yond the formal 5<r detection limit. In the more likely cases 
of Q decreasing as b^ 1 or 6~ 2 , b being the bias factor, the 
"clumps" will only show up in Planck maps as < 5<r fluctu- 
ations. Anyway, the Planck surveys will provide a rich cat- 
alogue of candidate proto-clusters at substantial redshifts 
(typically at z ~ 2-3), very important to investigate the 
formation of large scale structure and, particularly, to con- 
strain the evolution of the dark energy thought to control 
the dynamics of the present day universe. Detailed numeri- 
cal simulations carried out using the fast algorithm recently 
developed by Gonzalez-Nuevo et al. (2004) are fully consis- 
tent with the analytic results, although a full comparison 
would require an upgrade of the algorithm to include the 
effect of the evolving three-point correlation function. 

As the ratio of the beam-width to the clustering an- 
gular size decreases, the observed fluxes approach those of 
the brightest sources in the beam and the "clump" formal- 
ism no longer applies. However, simulations show that also 
in the case of the Herschel/SPIRE 500 pm survey the con- 
tribution of neighbours to the observed fluxes enhances the 
bright tail of the observed counts. Due to the extreme steep- 
ness of such tail, as predicted by the model of Granato et al. 
(2004), even a modest addition to the fluxes of the bright- 
est sources may lead to counts at flux densities ~ 100 mJy 
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several times higher than would be observed with a high 
resolution instrument. It should be noted that, in the case 
of strong clustering, the canonical 5<r detection limit (shown 
by the vertical dotted line in both panels of Fig. 2) does not 
frees the observed counts from the confusion bias. 
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